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We study the extended Kitaev-Heisenberg (EKH) quantum spin model by adding bond-dependent 
off-diagonal Heisenberg term into the original KH model, which was recently proposed to describe 
the honeycomb Iridates. A rigorous mathematical mapping of spin operators reveals the intrinsic 
symmetry of the model Hamiltonian. By employing an unbiased numerical entanglement renor¬ 
malization method based on tensor network ansatz, we obtain the global phase diagram containing 
eight distinct quantum phases. By using the dual mapping of spin operators, each of the individual 
magnetic phase in the global phase diagram can be clearly understood. At last, we show that a 
valence solid state emerges as the ground state in the quadro-critical region where multiple magnetic 
phases compete most intensively. 


PACS numbers: 75.10.Jm, 75.10.Nr, 75.40.Mg, 75.40.Cx 


The family of layered honeycomb iridium oxides 
Na 2 lr 03 and L^IrOs have attract ed g reat interest both 
experimentally and theoretically 11-17] due to its possible 

lid 0 . 


relevance of Kitaev spin liquid [18|- Earlier theoretical 
studies proposed that the exchange from direct d-orbital 
overlap of Ir ions may result in a Kitaev-Heisenberg (KH) 
model 0 - Extensive study of such model has shown a 
variety of fascinating phenomena, including an extended 
spin liquid phase and quantum phase transitions into 
several well-understood magnetic ground states. Exper¬ 
imentally, neutron scattering data revealed the presence 
of the zigzag phase in Na 2 lr 03 [I|-0|. However, it appears 
to be hard stabilize within the HK model; one may con¬ 
sider additional exchange paths [13] or further neighbor¬ 
ing hoppings 12 ]. Recently, a generic bond-dependent 
off-diagonal exchange term is included in the extended 
Kiteav-Heisenberg (EKH) model [16]. By using a com¬ 
bination of classical techniques and exact diagonaliza- 
tion, 120 ° and incommensurate spiral order, as well as 
extended regions of zigzag and stripy order are obtained. 
This EKH model can be regarded as the minimal model 
to describe the essential physics of honeycomb iridium ox¬ 
ides. However, the effect of bond-dependent off-diagonal 
exchange term is unclear and the obtained phase diagram 
is not clearly understood. 


In this letter, we perform a systematic study of the 
EKH quantum spin model on honeycomb lattice by us¬ 
ing the multiscale entanglement renormalization ansatz 
(MERA) [20] . We obtained a global phase diagram con¬ 
taining eight distinct quantum phases. A rigorous math¬ 
ematical mapping of spin operators reveals the intrinsic 
symmetry of the model Hamiltonian. In particular, by 
rotation of spin operators obeying certain rules, the addi¬ 
tional bond-dependent off-diagonal Heisenberg exchange 
term can be mapped to Heisenberg or Kitaev interac¬ 
tions and vice versa. So the original EKH model can be 
mapped to the model of same form but with modified 


coupling constants. By using this dual transformation, 
each of the individual magnetic phase in the global phase 
diagram can be clearly understood. For instance, some 
seemingly complicated magnetic orders can be mapped 
into simple antiferromagnetic (AFM) order or ferromag¬ 
netic (FM) order. Furthermore, both the finite size effect 
and the results for infinite system are discussed. At last, 
we show that a valence solid state emerges as the ground 
state in the quadro-critical region where multiple mag¬ 
netic phases compete most intensively. 

Model Hamiltonian .— Our model is the EKH model 
defined on a honeycomb lattice, for which the Hamilto¬ 
nian is 

H = Yi PS; • s i + KSJSJ + T(S?Sf + S?Sf )], (1) 

(i,3) 


where (i,j) denotes the nearest-neighbor, and terms with 
strength J, K and T represent Heisenberg, Kitaev and 
bond dependent off-diagonal exchange interactions, re¬ 
spectively. We use 7 to distinguish one spin direction 
on each bond for the Kitaev exchange, as shown in 
Fig. m We have parameterized the energy scales so that 
J 2 + K 2 + T 2 = R 2 , J = R cos 0 sin <f>,K = R cos 0 cos 0, 
T = RsinO. <f> controls the ratio between Kitaev and 
Heisenberg terms. 

It is well known that at the Kitaev limit J = T = 
0 ,K 7 ^ 0, the ground state of system is a topological 
spin liquid which has no long range order. Away from 
the limit, several long range magnetic orders may emerge 
and compete with each other. In addition to the AFM 
and FM order, the system also harbors stripy (ST) and 
zigzag (ZZ) phases, in region centered at two points K = 
±2 J. It can be shown that exactly at these two points, 
the original KH hamiltonian can be mapped into a pure 
Heisenberg model with the coupling constant J f = — J. 
The trick is to divide the system into four sublattices, 
each flip spin components in accordance to certain rules, 
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FIG. 1: (color online). Mapping the original EKH model to its 
counterpart through spin rotations. The upper panel stands 
for the original spin model, where spin directions 7 (a, /3) for 
Kitaev and bond-dependent (r) exchanges are denoted by 
different colors. By performing the spin rotation marked by 
numbers 0 to 5, following the rule listed in the table, we can 
map the original spin model to its counterpart. Different col¬ 
ors are painted for each site to denote rules of sign change, 
following which ZZ/SP phases can be mapped into AFM/FM. 
Shades in the background in the upper panel represent the 
MERA structure we adopted in numerical calculation. 

as denoted by different colors in FigJTJ In that sense, the 
ground states of KH model can be well understood in 
terms of the competitions among simple magnetic long- 
range orders. 

Mapping of Spin Operators. —In order to better under¬ 
stand the effect of T-term, we applied rotation of the spin 
operators S 7 -A S a (See the table in Fig.l) so that the T 
term can be mapped back to Kitaev exchange term. By 
doing so, the original hamiltonian is mapped to the model 
of same form but with modified coupling constants, 

H = • s 3 + K'SJS] + r '(sfsf + S?S?){ 2) 

<ij> 

where strength of each term is related to original param¬ 
eters as J r = —T, K' = {K — Y — J), and T' = —J. Notice 
that the Heisenberg and bond-dependent exchange cou¬ 
plings are swapped through the mapping. 

The mapping can be extremely helpful to understand 
the characteristics of ground state. For instance, we 
take the limit J = 0,T = K = 1. Under rotations of 
spin operators, the mapped Hamiltonian corresponds to 
a simple Heisenberg model with FM interaction H = 
— J2(i j) &iSj- Indeed, numerical calculation confirms 



FIG. 2 : Global phase diagram from MERA calculations of 
24-sites system. Longitude <f and latitude 6 of the world map 
parameterize the couplings of J, K and T terms. In the up¬ 
per panel, phases are described from the perspective of the 
original spin model, whereas the lower panel shows the phase 
map when the rotated spin model is considered. Regions of 
same type of magnetic order are shown in same color, la¬ 
beled in the bottom. Diamonds notates special points where 
the original model can be simplified to simple Heisenberg in¬ 
teractions. Red squares mark parameter sets for which we 
show spin configurations of ground states obtained from ER 
calculations. Phase boundaries marked by red wide lines cor¬ 
responds to first order phase transition. 

that its ground state is classical FM state with energy 
per-site exactly equal to —0.375. In Fig.[3j panel (A), we 
show details of the ground state, from perspectives before 
and after rotations of spins, for comparison. In terms of 
original representation, the spin configuration of ground 
state is rather complicated, and long range magnetic or¬ 
der can not be easily distinguished. On the other hand, 
in the mapped representation, the spin configuration of 
ground state is a simple ferromagnet with spins align¬ 
ing in same direction. Hence, the intrinsic symmetry 
of the model Hamiltonian may help us to better under¬ 
stand the seemingly complicated magnetic orders of the 
system induced by the additional off-diagonal exchange 
term. There exists several such parameter sets where 
the original model can be simplified to pure Heisenberg 
models (shown in Fig. [2]). These special cases offer some 
additional hints for solving EKH model, but it becomes 
rather complicated when various magnetic orders com¬ 
petes with each other, especially when the T term domi¬ 
nates. 

Global phase diagram .—We use an unbiased MERA al¬ 
gorithm to fully explore the phase diagram of EKH model 
on honeycomb lattice. In the upper panel of FigJIJ we 
show the ER tensor network we construct for the hon¬ 
eycomb lattice. This method is capable of capturing the 
Kitaev spin liquid and other phases where quantum fluc¬ 
tuations play an important role. The system sizes we 
calculated range from 24 to 96. 

The global phase diagram we numerically obtained is 
shown in Fig. [2j which is presented in format of a world 
map as functions of two parameters 6 and <f. The anti- 
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ferromegnet Kitaev point (spin liquid) is located exactly 
at the map center, and the original Kitaev Heisenberg 
model with T = 0 corresponds to the equator line. In 
Fig. [2j we present two phase maps and each of them be¬ 
longs to a different spin representation. The very same 
state show different magnetic long range orders without 
(upper panel) or with (lower panel) spin rotations. Here 
we denote these two different orders with labeling “X” 
and “F”. By performing the spin rotations, we notice 
that Heisenberg and T exchanges are swapped, namely, 
J = — T and vice versa. 

AFM/FM and 12(F phases —Due to the fact that the 
EKH model is mapped to itself, we expect the phase map 
to be self consistent by spin rotations. Such mapping cor¬ 
respondence is clearly observed for AFM and FM orders, 
which occupy a large portion of the phase map. It is im¬ 
portant to notice that one block of certain order (AFM 
for example) in the original map can be further sepa¬ 
rated into sections of different orders (120 and FM) in 
the phase map for rotated spins, vice versa. 

This feature seems to be contradictory, but careful 
analysis reveals that these two regions are indeed of dif¬ 
ferent magnetic orders. As a matter of fact, a first order 
phase transition is clearly shown when the phase bound¬ 
ary is approached and eventually crossed (see Figure in 
supplementary material). The physical reason can be 
understood as the following. For the FM order, spins 
(Si) are parallels for all sites. Naturally, we expect equal 
spin components for every site, namely Sf = Sf. On the 
other hand, the very phase is of the AFM order in term 
of original spins, which require \Sf\ = \Sf\. Notice that 

rotations of spins (Sf -A — Sf , for example) are involved 
when original hamiltonian is mapped to the new one. As 
a result, the only possible outcome for these two require¬ 
ments to be met at the same time is that all spins order 
along (1,1,1) direction, namely Sf = Sf = Sf. In an¬ 
other word, the bulk of FM/AF phase is of Ising type. 
This is to be expected, as introduction of T term breaks 
SU(2) symmetry of original Heisenberg magnets. 

Similar situation occurs for the 120°/AFM phase. The 
120° phase is named after the fact that all angles between 
next-nearest neighbor spins S are 2tt/3 in this order, an 
illustrated in Fig. 3. Such phase appear south of euqator, 
which represent a positive nearestjieighbor Heisenberg 
couplings J = —T, unfavoring the FM phase. Therefore, 
spins in original picture orders in the plane perpendicular 
to (1,1,1), as a compromise. (Notice, the system is still 
in the AFM phase with all nearest neighbors spins an¬ 
tiparallel.) It can be shown that whence Sf + Sf + Sf = 0 
(true for any vector in plane perpendicular to (1,1,1)), 
after rotations listed in tables, the 120° order emergies. 
In this case, the U(l) symmetry, instead of SU(2) sym¬ 
metry, is broken in the AFM order. 

Due to effect of T term, the orginal SU(2) symmetric 
Neel AFM order breaks into two different types: an Ising 
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FIG. 3: Spin configurations and nearest neighbors’ correla¬ 
tions (Si • S j) obtained from numerical calculation for 24-sites 
system. On each site, spin (Si) is plotted as combination of 
a circle (proportional to (Sf)) and a vector (proportional to 
( Sf y )). Spin correlations are notated by bonds connecting 
spin pairs, whose width are proportional to the strength of 
(Si • S j). Red/green color stands for FM/AF correlation, 
respectively. Panel A ~ F show results achieved in differ¬ 
ent phases, their corresponding parameter sets can be read 
from the marks on phase maps. The number of each nota¬ 
tion stands for the perspective from which the ground state is 
viewed: “1” notates the original spin picture, “2” stands for 
the rotated spin picture, and “3” means the result of zigzag 
mapping (for the rotated spin model), described in Fig. |T] 
Results for larger systems can be found in Sup. material. 


AFM north of the equator, and a U(l) AFM south. A 
first order phase transition shows up at the phase bound¬ 
ary. They are different physical orders with distinct bro¬ 
ken symmetry. Such phenomena is universal for all AFM 
and FM orders, in both original and rotated spin models. 
The finite size effect is less remarkable as shown in our 
numerical calculation. All magnetic orders remain the 
same, except a slightly shifts of phase boundaries [2l| . 

Zigzag/stripy phases and frustration —While switching 
longitude and latitude is helpful to understand the rela¬ 
tion between two phase maps, especially for AFM/FM 
orders, it does not work well when zigzag/tripy phases are 
considered. As a matter of fact, these phases are far more 
complicated than Neel-like orders discussed above. From 
numerical calculation, we find that spins are much more 
prone to tilting away than strictly order along one direc¬ 
tion. Meanwhile, various candidates of magnetic orders 













compete for the ground state, and energy gaps between 
them are relatively small. The reason for such different 
behaviors can be found in the mapping when rotations 
of spins are done. As mentioned before, such discrep¬ 
ancy may not affect AFM and FM orders much, (as we 
observed in our numerical calculation), but do change 
zigzag phase and stripy phase significantly. As a matter 
of fact, the new pattern of Kitaev and T interactions is 
frustrated in terms of zigzag/stripy order. 

For the original Kitaev Heisenberg model, a four- 
sublattice mapping between the ST/ZZ order and the 
AFM/FM order is well known and discussed [19]. Nat¬ 
urally, we expect same type of operation can be applied 
for the rotated EKH model. Unfortunately it is not pos¬ 
sible to obtain such a perfect mapping in this case. One 
can show that, when we flip spin operators one by one 
along an hexagon, according to arrangement of Kitaev 
interaction jij on each bond, we can not obtain self- 
consistency after a full circle. As a result, we have to 
sacrifice certairybonds in order to map the rotated spin 
model to its ZZ/ST counterpart. One way to make the 
compromise is shown in the lower panel of Fig. [lj We 
keep the sign-flipping rule along chains following certain 
direction, and reverse spin operators when jump between 
chains. Clearly there exists 3 degenerate ways (directions) 
for such sign change rules. Only one is plotted injffie fig¬ 
ure, and shades indicate the chosen direction of ZZ/ST 
chains. 

Such a compromise sacrifices certain Kitaev interac¬ 
tions intra-chain. Numerical study shows that such a 
mapping is energy-favorable for most zigzag region in 
the rotated spin model. Examples of such ZZ phases are 
shownMn Fig. [3] In the rotated spin model, we expect 
the ZZ phase to appear near the point K = — 2J, T = 0, 
which is translated to J = 0, if = — 1,T =1. Notice that 
the spin-spin correlation shows that intra-chains’ correla¬ 
tions are significantly larger than inter-chains’, which is 
consistent with the mapping we made. It is necessary to 
mention that due to this frustration, these states become 
less stable when we tune parameters away from the zigzag 
point, spin bend away and new magnetic order appears. 
One notice that for J < 0 and J > 0, the zigzag phase 
for rotated spin model also emerges as different magnetic 
order in original spin picture. The reason is very similar 
to the AF case explained above. Due to perturbation of 
T — — J term, system magnetize is different quadrant in 
these two cases. Details can be found in [2lj |. 

The mapping we choose is not necessarily the unique 
one, however, especially in larger systems. In small sys¬ 
tem with 24 and 36 sites, numerical calculations shows 
that the former mentioned states are the ground states 
spanning the zigzag region for the rotated spin model. 
However for larger systems, various magnetic order can¬ 
didates contradict to such mapping also appears, with 
reasonably good energy 0. With current calculation 



FIG. 4: Evolution of spin correlations (SC) (Si • S j) when 
system undergoes a transition from magnetic long range or¬ 
ders to the VBS phase. K and T are fixed to be 1, and J 
is tuned through the phase transition. Different colors repre¬ 
sents various distances between two spins i and j , from nearest 
neighbor (NN) to 5th nearest neighbor (5NN), the maximal 
distance in a periodic 24-site lattice. At — J/K — 0.95, a 
phase transition between FM and the VBS order occurs, and 
correlations for rotated spins (SC notated by suffix M) drops 
abruptly. When J/K goes over 1, the system enter the FM 
order, and SCs for original spins (suffix O) do not decay by 
increasing the separation. The inset shows the ratio between 
average spin correlation intra and inter a plaquette hexagon 
SC 1 /SC 2 ■ Results for various system sizes are plotted. 


capability, we are not able to probe the region thoroughly. 

Valence bond solid. — As expected, multiple magnetic 
long range orders compete most intensively close to two 
points — J = K = T = 1 and J = —K = —T = 1. 
These two special points corresponds to verging points 
of multiple phases, as shown in Fig. [2] It is not clear 
what is to be expected from analytical analysis, and our 
numerical calculation shows that small regions surround 
these multi-critical points harbor the VBS phase. A pe¬ 
riodic strong-weak shifting pattern for nearest neighbor 
spin correlation can be observed in Fig. [3j panel D2. 
More precisely speaking, correlations inside a plaquette 
(hexagon unitcell) is larger than intra-plaquette ones. As 
a result, we expect long range spin correlations to drop 
fast when the distance between two separated spins is 
increased. In contrast, for all magnetic orders discussed 
above, correlations decay much slower when spins are 
further separated. In Fig. [4j we show such results from 
numerical calculations and a transition between the mag¬ 
netically ordered phase (120/FM) and the VBS phase 
can be distinguished. 

To further clarify this, we define spin correlations in¬ 
tra/inter plaquette hexagon (P) as SC\ and SC 2 , where 
SCi = Ep,=p- Si • S, and SC 2 = £ P<#P . S* • S,, In 
subset of Fig. [4j we show that for various system sizes 
we calculated, the ratio between SC\ and SC 2 increase 
significantly when the VBS region is approached. It is 
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worth mentioning that the local spin expectation value 
(Si) vanishes in these critical regions (as shown in Fig. [3] 
for 24-sites system), a similar feature observed for the Ki- 
taev spin liquid. However, long range spin correlations 
reveal that they are intrinsically different. We have to 
clarify that our calculation is limited to relatively small 
system sizes, as a result, finite size effect boundary con¬ 
ditions may play important roles. We can not rule out 
the possibility of spin liquid in this critical region. 

Conclusions. —We numerically studied the global 
phase diagram of the EKH model containing eight dis¬ 
tinct quantum phases. A rigorous mathematical map¬ 
ping of spin operators reveals the intrinsic symmetry of 
the model Hamiltonian. By using the dual mapping of 
spin operators, each of the individual magnetic phase 
in the phase diagram can be clearly understood. Fi¬ 
nally, we showed that a valence solid state emerges as 
the ground state in the critical region where multiple 
magnetic phases compete most intensively. 
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